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We compute the continuum thermo-hydrodynamical limit of a new formulation of lattice kinetic 

^. ■ equations for thermal compressible flows, recently proposed in [Sbragaglia et al., J. Fluid Mech. 628 

299 (2009)]. We show that the hydrodynamical manifold is given by the correct compressible Fourier- 

Navier-Stokes equations for a perfect fluid. We validate the numerical algorithm by means of exact 

results for transition to convection in Rayleigh-Benard compressible systems and against direct 

f "•) ■ comparison with finite-difference schemes. The method is stable and reliable up to temperature 

rvj ' jumps between top and bottom walls of the order of 50% the averaged bulk temperature. We use 

this method to study Rayleigh- Taylor instability for compressible stratified flows and we determine 

the growth of the mixing layer at changing Atwood numbers up to At ~ 0.4. We highlight the role 

played by the adiabatic gradient in stopping the mixing layer growth in presence of high stratification 

r ~\ , and we quantify the asymmetric growth rate for spikes and bubbles for two dimensional Rayleigh- 

Taylor systems with resolution up to L x x L z — 1664 x 4400 and with Rayleigh numbers up to 

Ph ■ &~2x 10 10 . 

I. INTRODUCTION 

Lattice implementations of discrete- velocity kinetic models have gained considerable interest in the last decades, as 
efficient tools for the theoretical and computational investigation of the physics of complex flows [1-8]. An important 
^TJ ' class of discrete- velocity models for ideal fluid flows, the lattice Boltzmann models (LBM) [9-11] can be derived 
_f~. , from the continuum Boltzmann (BGK) equation [12], upon expansion in Hermite velocity space of the single particle 
distribution function, f(x,£,t), describing the probability of finding a molecule at space-time location (x, t) and with 
velocity £ [4, 13-15]. As a result, the corresponding lattice dynamics acquires a systematic justification in terms 
of an underlying continuum kinetic theory. The state-of-the-art is satisfactory concerning iso-thermal flows, even in 
presence of complex bulk physics (multi-phase, multi-components) [1, 2, 16] and/or with complex boundary conditions 
such as roughness, non- wetting walls and slip-length [6, 17-19]. 

The situation is much less satisfactory when temperature plays an active role in the flow evolution, due to complex 
compressible effects which are present even in ideal fluid/gas or to phase-change in multi-phase systems, or both. 
Only a few years ago, one could frankly admit that not a single known Lattice Boltzmann approach could handle, in a 
realistic way, thermal problems properly. The main difficulties being the development of subtle instabilities when the 
local velocity increases. In the last years, the situation has started to improve, with different attempts being made 
to describe active thermal modes within a fully discretized Boltzmann approach [15, 20-28]. These studies show that 
in order to recover the right continuum descriptions with the correct symmetries for the internal energy flux, one 
needs to enlarge the number of discrete speeds (a possible choice, for space filling schemes following a Gauss-Hcrmitc 
quadrature [15, 26], is 37 speeds in 2d [26, 29] and 107 speeds in 3d [30]), or to add ad-hoc counter-terms canceling 
spurious anisotropic operators [21, 22]. Otherwise, different hybrid attempts have been proposed, where temperature 
evolution is solved using finite difference methods [20] or with lattice schemes able to reproduce thermal Van der 
Waals fluids in the continuum limit [24]. Boundary conditions [23, 31] and stability issues [32] are also much more 
involved when thermal modes are present. It is fair to say that not a single model emerged as the optimal choice, and 
only a few explorative studies have been performed in order to check potentiality and limitations of each proposed 
solution. 

The aim of this paper is twofold. First, we intend to further discuss a recent formulation, proposed by some of us 
in [33], for a new way to incorporate the effects of external/internal forces in thermal LBM. We provide here the full 
explicit Chapman-Enskog expansion, whose results where only anticipated without proof in [33] , in order to show the 
convergence of the model to the Fourier-Navier- Stokes equations. We validate the method in a case where thermal 
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compressible effects play a major role, i.e. the transition to convection in a compressible Rayleigh-Benard system of 
height L z , with an imposed temperature jump, T u —Tj = AT. For such systems, it is possible to calculate the critical 
Raylcigh number analytically [34] at changing both the stratification parameter (also known as the scale height), 
Z = AT/T U , and the polytropic index, m = g/(R/3) — 1, where R is the gas constant, g the gravity acceleration and 
f3 = AT/ L z the temperature gradient. We show here that our LBM scheme is able to handle temperature jumps 
as high as AT/T U = 2 for both positive and negative values of the polytropic index (stable and unstable density 
stratification). Such systems are clearly very far from the classical Obcrbcck-Boussincsq approximation [35, 36]. 

Second, we study highly compressible Raylcigh- Taylor systems, for the initial configuration where two blobs of the 
same fluids are prepared with two different temperatures (hot, less dense, blob below, cold, denser, blob above). We 
show that the method is able to handle the highly non-trivial spatio-temporal evolution of the system even in the 
developing turbulent phase. In this case, we could push the numerics up to Atwood numbers At ~ 0.4. Maximum 
Rayleigh numbers achieved are Ra ~4x 10 10 for At = 0.05 and Ra ~2x 10 9 for At = 0.4. We present results on: 
(i) the growth of the mixing-layer at changing the compressibility degree, including the asymmetry in the quadratic 
growth of spikes and bubbles dynamics; (ii) a new effect of stratification which stops the mixing length growth when 
a critical width, L a d is reached. We interpret this as due to the existence of the adiabatic gradient: when the jumps 
between the two moving fronts leads to a temperature gradient, AT / L a d, of the order of the adiabatic gradient the 
dynamics stops and only thermal diffusive mixing may further acts. 

Technically speaking, the main novelty of the thermal-LBM formulation proposed in [33] relies on the fact that it 
is possible to incorporate the effects of an external and/or internal force (gravity and/or intermolecular potential) via 
a suitable shift of both momentum and temperature appearing in the local equilibrium distribution of the Boltzmann 
collision operator. Doing that, the systems acquires an elegant self-consistent formulation and a stable spatio-temporal 
evolution also in presence of compressible effects, as demonstrated by the examples anticipated before and detailed 
later. 

The paper is organized as follows. In Section II we briefly remind the details of the LBM formulation and we discuss 
the first result of this paper: the continuum thermo-hydrodynamical limit, given by the Fourier-Navier-Stokes eqs., as 
obtained from a rigorous Chapman-Enskog expansion of the discrete model. In Section III we show first the validation 
of the discretized algorithm by studying the transition to convection in compressible Rayleigh-Benard systems and 
comparing the results with exact analytical calculations, at changing the scale height and the polytropic index. In the 
same section we also present validation of the method against finite-difference methods, for the same set-up but after 
the transition, once convective rolls are present and stationary. In Section IV we show the investigations of another 
non-trivial compressible case: Rayleigh- Taylor system, for two different Atwood numbers At = 0.05 and At = 0.4. 
Conclusions and perspectives close the paper in Section V. 

II. THERMAL KINETIC MODEL AND CONTINUUM THEORY 

The main goal of this section is to show how to use a Thermal-LBM to discretize continuum thermal kinetic equations 
in presence of internal/external forces and how to extract via a suitable Chapman-Enskog multiscale expansion the 
relative hydrodynamical evolution, given in term of the forced Fourier-Navier-Stokes equations. The first issue was 
already discussed in [33]: here we briefly recall it and then discuss the second issue in details. 

A Thermal-Kinetic description of a compressible gas/fluid of variable density, p, local velocity w, internal energy, /C 
and subjected to a local body force density, g, is given in the continuum by the following set of equations (repeated 
indices are summed upon): 

dtp + dtipiii) = 

dt{pu k ) + di{Pik) = P9k (1) 

d t IC + \d % qi = pgtiit 

where Pik and qi are momentum and energy fluxes (still unknown at this level of description). 

In [33], it is shown that it is possible to recover exactly the above set of equations, starting from a continuum 
Boltzmann Equations and introducing a suitable shift of the velocity and temperature fields entering in the local 
equilibrium: f^ eq \£,; p,T,u) — > f^ eq \$,; p,T ', «). The new -shifted- Boltzmann formulation being: 

%+Z-Vf = -\{f- f (e?) ); (2) 



Where the shifted local velocity and temperature must take the following form: 

u = u + rg T = T- r 2 g 2 /D. (4) 

The lattice counterpart of the continuum description (2) can be obtained through the usual lattice Boltzmann dis- 
cretization: 

Mx + c t At,t+ At) - Mx,t) = _** (/,(*, i) - // e9) ) 

T 

where the equilibrium is expressed in terms of hydrodynamical fields on the lattice, /, (a;, p( L > ,vS L > ,T^ L '), and the 
subscript I runs over the discrete set of velocities, cj. The superscript L indicates that the macroscopic fields are now 
defined in terms of the lattice Boltzmann populations: 



,W = 

Dp{ L) T (L) =T \ Cl _ u iL)\* fl _ 



p (L) u (L) =Y.lCifl\ (5) 



In [33] it was shown that the lattice version of the shifted fields entering in the Boltzmann equilibrium (see Appendix 
A for its detailed form) is: 

uW = u w + Tg fw = rw + T{At D T)g2 + O(At) 2 . 

As it is known, lattice discretizations induce non trivial corrections terms in the macroscopic evolution of averaged 
hydrodynamical quantities. In particular both momentum and temperature must be renormalized by discretization 
effects in order to recover the correct thermal kinetic description (1) out of the discretized LBM variables. Density is 
left unchanged, p( H ' = p, while the first non trivial correction to momentum is given by the pre and post-collisional 
average [37, 38]: 

U {H) = U (L) + Ai g (6) 

and the first non-trivial, correction to the temperature field by [33] : 

T ( H ) = T (L) , ( At )V (7) 

4D ' w 

Using this renormalized hydrodynamical fields, one recover by a suitable Taylor expansions in At the thermo- 
hydrodynamical equations [33]: 

'd tP W+d l {puf ) )=Q 

dt(p (H H H) ) + di(PL H) )=P CH) 9k ( 8 ) 

d t KW + ld iq W=pW gi uf\ 

The above equations arc still unclosed. A closure ansatz to express the stress tensor, P> k , and the heat flux, q\ , 
in terms of lower order moments is needed. This ends our short review of the backup material. 

We proceed now with a systematic multi-scale closure of (8) in order to control the small wave-length limit where the 
full Fourier-Navier-Stokes equations emerge. The main added value with respect to previous similar calculations [40] 
is the explicit inclusion of the effects of the external force g in the Chapman-Enskog expansion. 
In order to perform the calculations, we need to introduce a hierarchy of temporal and spatial scales, via the intro- 
duction of a small parameter, e: 

d t -» ed u + ( 2 d 2 u di -> edi 

and the corresponding expansion for the Boltzmann distributions 

/ = /(0) + e/ (D + £ 2 / (2) + e 3 /( 3) + £ 4 /( 4) + 



together with a suitable rescaling of the forcing terms, g ~ 0(e) [37]. The various rescalings immediately reflect in 
the explicit expansion of the equilibrium distribution in terms of Hcrmite polynomials, n t : 



f° = ^£:U B) «, w 



where wi are suitable weights [27, 29]. The projections on the different Hermite polynomials, Oq , are explicitly given 
in Appendix A. 

After a long calculation, fully detailed in the Appendix, one shows that the leading long wavelength limit coincides 
with the continuum Fourier-Navier-Stokes equations of an ideal compressible gas given by: 



pdtuf* + puWdiuf) = P g 3 + %*$? (9) 



'd tP + d l (pu[ H) ) = o 

< pdtUj + pu\ OiU^ 
^e^ + puf > die & = <j\f ] d 3 uf J) + ^(fc^e^); 

with the ideal gas internal energy given by: e^' = S-T^- H \ The stress tensor is given by 



a\f ] = -pT^Sij + u(d iU ^ H) + djuf*) + <% U - f ) d k u 



(H) 
k ■ 



The shear and bulk viscosities are: 



•*»*-%): (^)-^h A 4 



and the thermal conductivity: 



T^p (r-^-). (10) 
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These are therefore the equations for a compressible gas with an ideal equation of state: 

P = P T {H) (11) 

and ideal specific heats: 

Cv = yj c p = - + 1 (12) 

It is not difficult to show that in the case the external forces are conservative, written in a potential form depending 
only on the density, one may easily incorporate these effects in the definition of an internal energy, opening the way 
to discuss also non-ideal equations of state [33] . 

III. TRANSITION TO CONVECTION IN RAYLEIGH-BENARD COMPRESSIBLE SYSTEMS 

A first non trivial application of the above algorithm can be found studying the behavior of Rayleigh-Benard cells 
both considering the effects of compressibility and stratification to the transition from diffusive to convective dynamics 
[34, 41, 42] or to the case of fully turbulent non-Obcrbcck-Boussinesq convection [36]. Here we concentrate on the 
first issue (see top panel of figure 1 for a schematic view), results on high Raylcigh turbulent convection will be 
published elsewhere. First, let us rewrite the set of equations (9) in a more transparent way, dropping for simplicity 
the superscript H in all variables and using the explicit expression of the internal energy in term of the temperature 
field: 

' 'D t p = -pdiUi 

pD t u t = -dip - pgS itZ + rjdjjUi + ( 1 - ^M rjdidjUj (13) 

^pc v D t T + pdiUi = kdaT + ri{d t Uj + djU z - j-5ijd k u k )diUj 

where we have introduced the material derivative, D t = dt + Ujdj, and we have assumed constant viscous and thermal 
conductivity coefficients [34, 43] . The equation of state is, p = pT, i.e. it is given in terms of quantities normalized such 



that the gas constant is R = 1. For a cell of height L z and with imposed bottom and top temperature, Td and T u , the 
hydrostatic equilibrium is easily found in terms of the temperature jump across the cell, [5 = {Td — T u )/L z = AT/L Z : 

'T {z) = {T d + T u )/2-pz 

p (z) = p(T (z)/f)™ (14) 

Po(z)=p(T (z)/f) m + 1 

where the two integration constants must satisfy, p — pT, with T a reference temperature, T = (T„ + Td)/2. In 
(14) we have introduced also the polytropic index: m = g/j3 — 1. At changing the polytropic index, one changes 
the hydrostatic profiles of density and pressure. In order to be unstable, the profile must obviously verify, /3 > 
(if g > 0, as assumed here) and therefore the interesting polytropic interval is limited to m > — 1. Furthermore, 
unstable fluctuations may develop only if the hydrostatic temperature gradient, (3 is larger than the adiabatic 
gradient, f3 ac [ = g/c p , i.e. only when the adiabatic transformation of a hot/cold spot of fluid moving up/down 
induces a temperature variation that does not exceed the hydrostatic change [44] . This limits the interesting interval 
excursion of the polytropic index from above, m < c p — 1, which in our units, for an ideal gas in 2d, means m < 1. 
The limitation from above is a typical important example induced by compressibility/stratification, i.e. by the fact 
that a cold/hot fluid spots may contracts or expand during their spatio-temporal evolution. Stratification can be 
also measured by the scale height, i.e. a typical length scale, Lh, built in terms of mean hydrostatic quantities. In 
our case, the most natural way to define it is by using the temperature profile: Lh = (T U /AT)L Z = L z /Z. Where 
we used the dimensionless parameter, Z = AT/T U which is a direct measurement of the stratification effects: for 
Z ^> 1, the cell height L z is much larger than the typical stratification length, i.e. the fluid is highly stratified. 
On the other hand, the limit Z — > corresponds to the so-called Oberbeck-Boussincsq approximation, where 
both stratification and compressibility are vanishingly small. The latter is, by far, the most studied convection 
configuration, even though some important applications for astrophysics [45, 46] and recently also for laboratory 
set-up [47-49] cannot neglect compressible modes. It is possible to show [35] that in the Boussinesq approximation, 
the dependency from the polytropic index disappear (as it must obviously do) while it remains a possible effect 
induced by the adiabatic gradient (usually small on laboratory experiments, but not necessarily on atmospheric scales). 

We use this complex set-up to benchmark the thcrmal-LBM algorithm proposed, and probe its robustness at 
changing compressibility. This can be done directly against exact results on the emergence of convective instability 
in the system. It is possible to calculate, in a closed form, the stability problem of the linearized system around the 
hydrostatic solution (14), for both slip or no-slip velocity boundary conditions and for any polytropic index [34]: these 
are just suitable extensions of the well known Raylcigh calculation made for the incompressible case [50]. 
Stratification makes the problem non-homogeneous (in the vertical direction) and therefore it is not possible to define 
in a unique way the Rayleigh number. Anyhow, it turns out that it is possible to introduce a height-dependent 
Rayleigh number which rules the linearized problem: 

(g/T (z))LJ(f3 - p ad ) 
Ka[Z) ~ (k/p (z)c p )(W P o(z)) ' (15j 

and one can express the whole bifurcation diagram in terms the value of the Rayleigh number at a given height, 
say the middle of the cell z = L z /2 for example: Ra = Ra(L z /2). Different works have been devoted to the 
calculations of the critical Ra c at changing the polytropic index, the scale height, Z and the boundary conditions at 
the top/bottom plates [34, 51, 52]. A result of the stability calculation predicts that there exists a critical Rayleigh 
number which depends only on the polytropic index, m, on the stratification parameter, Z , and on the wavelength, 
a, of the perturbation, Ra c (m, Z,a). The hydrostatic solution will therefore become unstable under perturbation of 
a wavelength corresponding to the minimum possible critical Raylcigh number. Compressibility and stratification 
may have different effects, either stabilizing or destabilizing the systems, depending on the hydrostatic underlying 
equilibrium. For example, if the hydrostatic profile has an unstable density profile, m < 0, one gets that the 
critical Raylcigh decreases at increasing temperature jumps. The opposite happens when density is stably stratified, 
m > 0. From the definition of Raylcigh given in (15), it is easy to realize the importance of the adiabatic gradient, 
Pad = g/cp, i.e. if /? < I3 a( i, the control parameter is always negative and the system will always be linearly stable. 
In figure 2 we show the result of a numerical search of the critical Rayleigh number (i.e. the onset of the transition 
to convection) using our LBM algorithm, obtained by exploring the long time behavior of the system, prepared with 
a small perturbation to its hydrostatic equilibrium, and monitoring the successive temporal growth/ decline of the 
total kinetic energy (example in the inset). The LBM has been applied by imposing no-slip impenetrable boundary 
conditions for the velocity field at top/bottom walls, u z (z = ±L z /2) = 0; u x (z = ±L z /2) = 0; and with an imposed 
constant temperature jump, T(z = —L z /2) = Td', T(z = L z /2) = T u . Lateral boundaries are fully periodic. Technical 



details on the way to implement the given boundary conditions in the LBM algorithm are given in Appendix B. In 
the same figure we also report the critical Rayleigh numbers obtained from the LBM exploration, compared with the 
exact analytical results obtained by solving numerically the eigenvalue problem for the linearized equations as given 
in [34]. As one can see, the agreement is good, even for large temperature jumps, up to Z ~ 2. Larger values of Z are 
difficult to reach, because of limitations imposed by numerical stability of the boundary conditions and by the growth 
of unstable compressible modes in the system. In order to overcome such limitation one should probably extend 
the Hcrmite projections to higher and higher orders [32]. The main error source in the determination of the critical 
Rayleigh number out of our LBM method stems from the presence of spurious, small, departure from the exact linear 
profile in the mean temperature close to the boundary walls. This departure goes together with the existence of small 
spurious transverse velocity for two-three grid layers close to the wall and are due to the existence of discrete velocities 
which connects up to three layers in the lattice inducing non-local boundary conditions effects (see appendix A and 
B for details). Such effects can be annoying for the investigation of highly turbulent regimes, where the boundary 
layer dynamics becomes crucial to drive the correct thermal exchange with the bulk [53]. This shortcoming can be 
strongly reduced by moving from LBM algorithms using exact streaming (as done here) to LBM based on finite- 
volume schemes [54]. Details in this direction will be published elsewhere. The small spurious oscillations close to 
the boundaries does not prevent to get a very good quantitative validation of the algorithms also when large scale 
convectivc rolls are present. For example in figure 3 we make a one-to-one comparison of the LBM numerics with 
a numerical study using finite-difference scheme for incompressible Rayleigh-Benard systems [55, 56]. Again, the 
stationary profiles are perfectly superposing, as shown for both temperature and velocity in figure 3. This ends our 
validation section. In the next section we apply the new algorithm to study compressible dynamics, as it is the case 
of Rayleigh- Taylor instabilities in thermal stratified flows. In the latter case, the small spurious oscillations close to 
the walls are obviously completely unimportant, being the bulk the only physically interesting region. 

IV. RAYLEIGH- TAYLOR SYSTEMS 

Superposition of a heavy fluid above a lighter one in a constant acceleration field depicts a hydro-dynamic unstable 
configuration called the Rayleigh- Taylor (RT) instability [50] with applications on different fields going from inertial- 
confinement fusion [57] to supernovae explosions [58] and many others [59]. Although this instability was studied 
for decades it is still an open problem in several aspects [60]. In particular, it is crucial to control the initial 
and late evolution of the mixing layer between the two misciblc fluids; the small-scale turbulent fluctuations, their 
anisotropic/isotropic ratio; their dependency on the initial perturbation spectrum or on the physical dimensions of 
the embedding space [61, 62]. In many cases, especially concerning astrophysical and nuclear applications, the two 
fluids evolve with strong compressible and/or stratification effects, a situation which is difficult to investigate either 
theoretically or numerically. Here, we concentrate on the large scale properties of the mixing layer, studying a slightly 
different RT system than what usually found in the literature: the spatio temporal evolution of a single component 
fluid when initially prepared on the hydrostatic unstable equilibrium, i.e. with a cold uniform region in the top half 
and a hot uniform region on the bottom half (see bottom panel of figure 1). For the sake of simplicity we limit 
the investigation to the 2d case. While small-scales fluctuations may be strongly different in 2d or 3d geometries, 
the large scale mixing layer growth is not supposed to change its qualitative evolution [63, 64]. A grey-scale coded 
snapshot of a typical RT run is shown in figure 4 showing all the complexity of the phenomena. Let us start to define 
precisely the initial set-up. We prepare a single component compressible flow in a 2d tank of size, L x x L Zl with 
adiabatic and no-slip boundary conditions on the top and bottom walls, and with periodic boundary conditions on 
the vertical boundaries. For convenience we define the initial interface to be at height z = 0, the box extending up to 
z = L z /2 above and z = —L z /2 below it (see figure 1). In the two half volumes we then fix two different homogeneous 
temperature, with the corresponding hydrostatic density profiles, po, verifying [65]: 

dzPo(z) = -gpo(z). (16) 

Considering that in each half we have po(z) = T po(z) , with T fixed, the solution has an exponentially decaying behavior 
in the two half volumes, each one driven by its own temperature value. The initial hydrostatic unstable configuration 
is therefore given by: 




T u ; p (z) = p u cxp(-g(z - z c )/T u ); z > 

T d ; p (z) = p b exp(-g(z - z c )/T d ); z < 0. 



(17) 



To be at equilibrium, we require to have the same pressure at the interface, z — z c — 0; which translates in a simple 
condition on the prefactor of the above expressions: 

PuT u = p b T d . (18) 



Because T u < T<j, we have at the interface p u > pb- As far as we know, there are no exhaustive detailed calculations 
of the stability problem for such configuration, even though not too different from the usual RT compressible case 
[50, 66, 67]. As said, this is not the common way to study RT systems, which is usually meant as the superposition 
of two different miscible fluids, isothermal, with different densities [50, 61, 66, 68]. As far as compressible effects are 
small, one may safely neglect pressure fluctuations and write - for the case of an ideal gas: 

and the two RT experiments are then strictly equivalent. Moreover, in the latter case, if one may neglect the 
dependency of viscosity and thermal diffusivity from temperature, the final evolution is indistinguishably from the 
evolution of the temperature in the Boussincsq approximation [62, 63]. Here we will study both the case of small 
compressibility and small stratification, where pressure is always close to its hydrostatic value, p ~ po, and the case 
when compressibility becomes dynamically relevant, changing the global large scale evolution of the mixing layer. 

A. RT instability in thermally active flows: the role of the adiabatic gradient 

The main novelty in the set up here investigated is due to the presence of new effects induced by the adiabatic 
gradient, which in our case can be written as in the previous Section (3 ac i = g/c p . In order to understand the main 
physical point it is useful to think at the RT mixing layer as equivalent to a (developing) Rayleigh-Benard system with 
an imposed mean temperature gradient [69, 70]. Let us denote with L m i(t) the typical width of the RT mixing layer 
at a given time as measured for example from the distance between the two elevations where the mean temperature 
profile is 1% lower or higher then the bottom and top, respectively, unmixed temperature values, L m i — z u — Zd, where 
(T(x, z u )) x = 1 .01T„ and (T(x, Zd)) x — 0.99Td. It is well known that the temperature tends to develop a linear profile 
inside the mixing region, the resulting instantaneous temperature gradient is then given by /?(£) = (Td — T u )/L m i(t), 
and it decreases in time inversely to the growth of the mixing length. As a result, soon or later (if the box is tall 
enough) the instantaneous temperature gradient will become of the same order of the adiabatic gradient, j3(t) ~ (3 a d 
and the growth of the mixing length will stop. One can define an instantaneous Rayleigh number, driving the physics 
inside the mixing layer, estimated as in Section III: 

Rait) = (gz^imizM, (20) 

(k/ P oc p ){is/po) 

where (•) indicates quantities evaluated at the middle layer. It is clear that for small times, /3(£) <C (3 a d, the effective 
instantaneous Rayleigh number is high: the system is unstable, and the mixing length grows. On the other hand, 
as time elapses, the vertical mean temperature gradient decreases, until a point when, j3(t) ~ /3 a d, the instantaneous 
effective Rayleigh number becomes Ra(t) ~ 0(1) and the system tends to be stabilized. We can then identify an 
adiabatic length: 

Lad = (T d - T u )/fi ad = c p AT/g 

which determines the maximum length achievable by the mixing layer, in our configuration. Let us notice that in 
absence of the adiabatic gradient, the Rayleigh number would continue to grow indefinitely, being proportional to the 
third power of L m i(t), as it is the case for usual RT systems. If the profile coinciding with the adiabatic gradient is 
going to be fully stable depends on the top/bottom boundary conditions imposed on the whole spatial domain. In 
any case, when temperature matches the adiabatic profile, the system strongly feel it, showing a sudden slowing down 
of the mixing layer growth. To our knowledge, this effect has never been predicted before, within this framework. We 
show in figure 5 the evolution of temperature profiles when adiabatic effects are important. It is clear how the mixing 
layer growth is strongly slowed down when L m i(t) ~ L a d] afterward only very slow relaxation process happens further, 
mainly at the border between the edge of the mixing layer and the fluids region with homogeneous temperature. 
A possible way to estimate quantitatively when and how the adiabatic gradient starts to play a role in the growth of 
the mixing length is to use a simple phenomcnological closure for large scale quantities in the system. We start from 
the self-similar scaling predicted by [71, 72] for the homogeneous not stratified growth: 

(L ml (t)f = 4a< L > g At L ml (t) (21) 

which has a unique solution (beside the trivial one, L„a = 0) in terms of the initial value, L m i(to)'- 



L m l(t) = L ml (t ) + 2^L ml (t Q )a^Atg(t - to) + a (L) Atg(t- t Q ) 2 . (22) 



Eq. (21) offers the advantage to be local in time, i.e. one may extract the value of a^ by a simple evaluation of the 
plateau in the ratio (L m i) /L m i, time by time. In order to minimally modify the above expression considering the 
saturation effects induced by stratification, we propose to use: 



(L ml (t)f = Aa^gAtL ml (t)i, (*^® 



(23) 

where ip = ip(x) must be a function fulfilling the condition ip — >• 1 as x — > (that is for L a d — > oo), in order to 
recover the equation (21) for the not stratified case when the adiabatic gradient goes to zero. We further add the 
requirement of reaching the adiabatic profile with zero velocity and acceleration, enforcing a strict irreversible growth, 
i.e. L m i > 0, as it must be for the case of miscible fluids. Under these assumptions, it can be shown that the simplest 
form for the function ip is: 



L a d 



2Lgd — L 

Lad 



(24) 



where the prefactor C must be set equal to l/(e — 2) to comply with the prescribed boundary conditions. Equation 
(23) must be considered as a zero-th order phenomenological way to take into account of the adiabatic gradient in 
the mixing layer evolution. 

We integrated numerically eq. (23) testing the result in figure 6 where we show that it is possible to fit the 
global evolution of the mixing length L m i(t), by using reasonable [60] values of a^ L \ for all times, including the long 
time behavior where L m i(t) ~ L a d- In the same figure, we also show the behaviour of the time-dependent effective 
Rayleigh number (20), estimated using the instantaneous mixing length, L m i(t). As one can see, after the initial 
monotonic growth of the turbulent intensity, there appear a sudden slowing down, as identified by a strong reduction 
in the effective Rayleigh number. We can therefore safely assume that the solution of our equation (23) is a good 
generalization of (22) including also the adiabatic gradients effects. 

B. Compressible effects and mixing layer growth 

As shown in the previous section, effects induced by the adiabatic gradient start to appear when the mixing length 
becomes of the order of the adiabatic length L m i(i) ~ L a d- It is nevertheless possible to study the limit L{t) <C L a d 
but still observing important effects due to compressibility. Indeed, compressibility due to stratification is controlled 
by the Atwood number. From the expression of the instantaneous Rayleigh number (20) one may compute the 
typical length scale at which turbulence will be maximal, i.e. the largest extension of the mixing layer up to which 
the Rayleigh number is still growing, before decreasing because of the adiabatic gradient. This is just given by the 
maximum of Ra(t) as a function of time, which is reached at a characteristic time, t* such that: 

L ml {t*) = \L ad =^^. (25) 

It is also possible to estimate the typical Mach number reached at the maximal turbulent intensity, considering that 
hydrodynamical velocities can be estimated as, V ma x ~ d/dtL m i{t*) = 2a^ At g t* and that the minimal sound speed 
is given, in our units, by v s — \/Td, we get for the Mach number at the maximal turbulent intensity: Ma ~ AtyoM^Cp 
where we have used (22) to estimate t* at a given L m i{t*). As a result, dynamical compressibility is only driven by 
the Atwood number -at fixed c p . Using the typical values of a ( - L - ) ~ 5 10~ 2 , as reported in the literature [60], and 
plugging the correct prefactor, we estimate Ma ~ 0.4, for the largest Atwood we could achieve, At ~ 0.4. 
It is well known that compressibility effects break the up/down symmetry in the mixing layer propagation [71, 72], 
downwards spikes (cold fluid blobs) move faster than upwards bubbles (hot fluid blobs). Such effect is completely 
missing in Boussinesq approximation where there is a perfect up/down symmetry, by definition. 

Neglecting slowing down effects induced by the adiabatic gradients, i.e. limiting the study of the mixing layer growth 
up to L m i(t) <C L a d, we may investigate the symmetry breaking in our set up at changing the Atwood number. To 
give an idea of the effects of compressibility, we show in figure 7 a few instantaneous mean profile of temperature, 
density and pressure for the two Atwood numbers here investigated. From the density and temperature profiles it is 
easy detectable, already by naked eyes, the asymmetry present for the high Atwood case At = 0.4 in the growth of the 
mixing layer, with the colder and denser front moving faster. Also, the appearance of non-trivial fluctuations in the 
pressure around the hydrostatic profile, for the case at At = 0.4, are the clear evidence of compressible effects at play. 
Both the asymmetry and the pressure fluctuations arc completely absent for the case at small Atwood (left panels 
of figure 7, an evidence of Boussinesq- like thermal fluctuations). All numerical experiments have been performed by 



preparing the initial configuration in its hydrostatic equilibrium (17) plus a smooth interpolation between the two 
half volumes in order to have a finite width of the initial interface. The initial temperature profile is therefore chosen 
to be: 



To(z) 



T u + T d 



T u -T d ,/(-?- z c ) 
tanh ' 



where with w we define the initial width of the interface and z c its unperturbed height (z c = in our frame of 
reference). Initial density po(z) and pressure po(z) are then fixed by solving the hydrostatic equation (16) in order to 
get the hydrostatic solution corresponding to the smoothed temperature profile. 

To destabilize the initial configuration, we follow [73] and shift randomly the center of the interface by adding 
horizontal perturbation at different wavelengths in the range k € [k m i n : k max \: 



k=k ma 
k=k mi , 



cos(2ir k x/L x + <pk) 



(26) 



where <f>k are random phases and N = \Jk max — fc m m, in order to have a total amplitude for the initial width almost 
independent on the number of modes. We have tried different ranges of wavelengths, without observing quantitative 
differences in the large time growth of the mixing layer. The ratio W ~ e/w gives the "wiggling" of the interface, i.e. 
how much the perturbation of the interface position is important with respect to the interface width. 
Below, we present results in different geometry, up to a resolution of L x x L z — 1664 x 4400 with different choices of 
W . For each parameters set we made typically O(50) separate RT evolution, starting from different random phases 
initial configurations. 



In the sequel, we show a summary of the results from two typical numerical series of runs, one with At = 0.05 (small 
compressibility) and a second one with At = 0.4 (large compressibility). It is useful to adopt a different definition for 
the mixing length in terms of a bulk mixing percentage, introducing the characteristic function (tent-map): 



xK] = 2£; < £ < 1/2 
xK] = 2(i-0; i/2<£<i 



(27) 



and defining the mixing length as [72] 



m 'hi 



dxdzx 



T(x, z) - T u 
T d -T u 



(28) 



It is easy to realize that if the temperature is fully homogenized in the fluid, T(x, z) = (T u + T d )/2, then the mixing 
length coincides with the full vertical extension of the box: H = L z ; if we have two perfectly separated hot and cold 
regions we have H = 0. In the intermediate situation when we have a mean linear temperature profile for z £ [z dl z u ], 
between two unmixed regions (T = T u if z > z u and T = T d if z < z d ) the mixing length estimated by (27) is exactly 
given by half of the linear region, H — (z u — z d )/2.. The definition of the mixing length (27) must be preferred 
with respect to more common definition of L m i based on thresholds on the linear profile, as adopted in the previous 
section. The former, being based on a bulk measure is not affected too much on the highly fluctuating properties of the 
interface between mixed and unmixed fluids. This is particularly important in 2d, where the averaged profile, being 
a one-dimensional cut, may fluctuate a lot (see also figure 7). Anyhow, in the case of a perfectly linear temperature 
profile the two lengths are obviously related by the relation H = l/2SL m i, where S is the percentage threshold used 
to identify the mixing front (in the previous section 6 = 0.99). 

Moreover, because here we want to distinguish the downward growth of the front due to cold spikes from the upward 
growth of bubbles, we introduce two different integral mixing lengths: 



z X 



H.(t) = ^- f dxdz (e(f - 

H b {t) = j- J dxdz fei«-^] 



T(x,z) -T u 

T d — T u 
'T(x,z) -T u 



T d — T u 



where of course, H (t) = H s (t) + Hb(t). Clearly, the a^ value ruling the long term quadratic growth of the integral 
mixing H is not necessarily the same of L m i. Typically one expects the same relation cS H > = Q.5 5a^ L > valid for the 
definition of the two mixing length, at least for times long enough. 
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As one can see in figure 8 there is a wide scattering of the mixing length evolution from run to run, where the only 
differences between them is the realization of the initial random phases. Due to the intense local temperature and 
density fluctuations, averaging over horizontal direction is not very efficient to smooth down statistical fluctuations, 
and one observes high variations from sample to sample: many realizations are needed to extract stable quantitative 
results on the long time evolution. In order to have an insight on the typical fluctuations we decided to analyze run 
by run and following two fitting procedures. First, we start from the equivalent of (22), written for bubbles and spikes 
separately: 

(H b (t - t Q ) = H b (0) + H b (0) t + a[ H) Atgt 2 
1 H s (t - t ) = H 3 (0) + H s (0) t + a { s H) Atgt 2 



with Hh,s(0) = 2 y H btS (0) a b s Atg, where to must be understood as the time when the initial perturbation is fully 
entered in its non-linear regime. In other words, to must be larger than the typical characteristic time of the slowest 
unstable mode. It can be estimated from linear stability analysis as t ~ y / L x /(2TT g At). A brute force way to extract 

the growth rate is to evaluate the ratio a s b = lim^oo H sb (t)/t 2 . Even, neglecting possible contamination due to 
stratification, this is of course valid, only asymptotically, when both dependencies on the initial time to and on the 
initial mixing length H Sjb (to) become negligible. As a matter of fact, taking into account also the maximum time 
achievable due to numerical limitations, it is very difficult to extract stable statistical results on the a^ fluctuations 
starting from the brute force analysis of (29). For instance, we found that a parabolic fit to our data, taking a s b free 
is very sensitive to the initial time to and/or the initial distance H btS (to), without allowing for a systematic assessment 
of the asymptotic behaviour. To give an idea of the importance of the initial condition versus statistical fluctuations, 
we show in the bottom panel of figure 9 the results of the asymptotic ratio H sb (t)/t 2 for two different scries of runs 
with different initial conditions. As one can see, even if asymptotically there is a clear tendency to forget the initial 
separation, in agreement with (29) there is not a well developed plateau, up to the time achievable in out numerics, 
indicating the existence of important sub-leading effects. The existence of such terms is highlighted in the inset of 
the same panel, where a log-log plot clearly shows the lack of a plateau even for large times. 

Another alternative, and more robust way, to extract a^ relies on the differential equivalent of (29) given by (21) or 
(23) when stratification becomes important. Using (21), one may directly assess the non-linear growth rate, without 
spurious contamination from initial conditions. 

In the upper panel of figure 9 we show the same data plotted in the lower panel but for the ratio 

a { $ = {H s , b {t)) 2 /{AgAtH s , b {t)l (30) 

i.e. we address time-by-time the part depending on asymptotic growth rate only. It is evident the net improvement in 
both the extension of the range where aS H ' coefficients are constant and the clear disentanglement of effects coming 

from the initial conditions. Out of the data for I H s ^ b (t) ) /4 (g At H Sib (t)) we may estimate the statistical fluctuations 

of q^ b , by making a fit to a constant in a given time windows. In figure 10 we plot the results of fitting the evolution 
(30) independently for bubbles or spikes (upward or downward fronts). From this we learn a few interesting facts: 
(i) at small Atwood (upper panel) bubbles and spikes travels almost with the same statistics, even though a small 
asymmetry can be observed in the shape of the whole histogram. The asymmetry is so small, that if averaged 
quantities are measured, the differences between them falls within error bars; (ii) there are not important effects form 
initial conditions -compare the two upper panels obtained with two different classes of initial conditions-; at least 
when data are fitted using (30), confirming that the observed spatio-temporal evolutions is dominated by strongly 
non-linear fully developed dynamic; (iii) at large Atwood (lower panel) the asymmetry becomes evident, spikes are 

( }-f\ ( M\ 

systematically faster then bubbles, the two evolutions gives different mean vales for a s and a b parameter. Our 
measure of the average global growth rate a^ H ' , can be estimated by summing up the growth rate in the two half 
cells: a^ = a s + a b ~ 0.02 is agreement with values typically found in literature [60, 71, 72]. For instance, 
in [60] a detailed overview of numerical results gives for the growth rate of bubbles, measured on the 99% width, 

a b ~ 0.025 ± 0.003, in agreement with a b = 0.0095 ± 0.002 we found for our integral growth rate (see caption 
of figure 10) taking into account that by definition one expects a factor two between the measurement made on the 
integral quantity, o^H), and the measurement made on the 99% level set, eS L \ 

The last issue we want to discuss concerns with homogenization inside the mixing layer. It is easy to show that in 
the Boussinesq approximation for a convective stationary cell with a mean linear temperature profile, all deviations 
from the mean profiles are homogeneous. The case of RT evolutions investigated here is slightly different. First, 
whenever stratification is important, there is no reason to expect exactly homogenization inside the mixing length. 
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Second, and more importantly, homogeneity must be expected only well inside the mixing layer, far from the up and 
downside fronts, where clearly strong non-homogeneous effects for both mean and fluctuating quantities must appear. 
It is interesting therefore to test, how homogeneous the statistics is, also to quantify the degree of mixing. In order 
to do that, we introduce the p-th order moments of temperature fluctuations: 

QV>(z,t) = ({T{x,z)-{T{x,z)) x Y) x . (31) 

In figure 11 we show the root mean square fluctuations around the vertical mean temperature profile, Q^ 2 '(z,t), 
(bottom panel) and the flatness, F(z,t) = Q^'(z,t)/(Q^ 2 '(z,t)) 2 , i.e. the ratio between fourth and squared second 
order moments of fluctuating quantities (top panel). As one can see, the root mean square fluctuations tend -very 
slowly- to develop a flatter and flatter plateau inside the mixing region, demonstrating that if the mixing layer is 
wide enough, there will be a larger and larger region where statistics is pretty homogeneous. On the other hand, if 
we plot the Flatness as a function of a normalized mixing length width, it converges towards a self-similar profile, 
for any time, where the effects coming from the two boundaries of the mixing regions are felt inside the whole layer, 
without showing any trend towards homogenization. This second finding is a clear indication that if normalized with 
the total mixing length extension, the region where the statistics may be considered homogeneous does not increase 
with time. 

V. CONCLUSIONS AND PERSPECTIVES 

We have explicitly computed the continuum thermo-hydrodynamical limit of a new formulation of Lattice Kinetic 
equations for thermal compressible flows, recently proposed in [33] We have shown that the hydrodynamical manifold 
is given by the correct compressible Fourier-Navier-Stokes equations for a perfect fluid. We have validated the calcu- 
lations against exact results for transition to convection in Raylcigh-Bcnard compressible systems and against direct 
comparison with finite-difference methods. The method is stable and quantitatively reliable up to temperature jumps 
between top and bottom walls (stratification) of the order of AT jT u ~ 2. We have also applied the method to study 
Rayleigh- Taylor instability for compressible stratified flows and we determined the growth of the asymmetric mixing 
layer at changing Atwood numbers up to At ~ 0.4 and to Rayleigh Ra ~2x 10 10 . We determined the distribution 
of the growth rate for bubbles and spikes, at changing At and we discuss its dependence on the initial perturbation. 
We also discussed the importance of the adiabatic gradient for the growth of the RT mixing layer in strongly stratified 
systems. In the latter case, we showed the existence of a maximal width, the adiabatic length, L ac [, for the mixing 
region. The high flexibility -and locality- of LB algorithm makes them the ideal playground where to push the resolu- 
tion, having perfectly scalable performances as a function of the number of processors in the parallel architecture. In 
particular, it is simple to extend such algorithm to deal with fully 3d systems for ideal, non-ideal and/or even immis- 
cible two fluids systems. High resolution studies of Rayleigh- Taylor systems meant to investigate short wavelengths 
scaling properties of velocity, density and temperature fields for high Rayleigh, with and without surface tension [39] , 
and using a highly optimized LBM algorithm for the Cell Broadband Engine [75] are under current investigation 
and will be reported elsewhere [76]. The thermal LBM here proposed still suffers of small spurious oscillations of 
temperature and perpendicular velocity close to the solid boundaries, making it still not appropriate to study high 
Rayleigh numbers stationary convection. A possible way to overcome this difficulty consists in abandoning numerical 
schemes based on exact streaming and to develop the proposed thermal LBM on a finite volume scheme. Results in 
this direction are out of the scope of this paper and will be the subject of a forthcoming publications. 
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VI. APPENDIX A 

In this appendix we detail the steps of the Chapman Enskog expansion leading to the thermohydrodynamical 
equations under the effect of general forcing term pg. Similar analysis (without the effect of the forcing) can be found 
in [40]. We start from the shifted equilibrium formulation 

fi(x + ciAt, t + At)- f t (x, t) = ~ (fi(x, t)-fi) (32) 
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where, for the sake of simplicity, in the notation of this appendix we have renamed the equilibrium distribution 



function with shifted fields, /, 



(eq) 



fr 



fi = fi(pM L) +X,T {L) +X) 



and where x an( i X are general momentum and temperature shifts for the equilibrium distribution with u^ L \T^ L ' the 
lattice velocity and temperature hereafter denoted simply with u and T. Central to our analysis is the expansion of 
the equilibrium distribution in Hermite polynomials [15, 26, 40] 



fl 



wi 



^^4 n) ( P ,« + £X,T + £ 2 A)^ 



(») 



with wi suitable weights whose values are reported in [26, 29] for the D2Q37 model here used (sec also figure 12). For 
the purposes of our investigation a fourth order approximation proves to be enough to recover the correct equations 
with the right isotropic properties for all hydrodynamical fields and tensors up to the eighth order [29] . The Hcrmitc 
polynomials are given by the following relations: 



H 



(o) 



H 



(i) 



cv, 



H 



(2) 



(33) 



uf ] =c 3 -8 Cl - 



n { , 4) =cf- 8c? 



ss 



(34) 



(«) 



and the projection coefficients a\ by 



,(°) 



(i) 
Qq = pu 



£PX 



,(2) 



,(3) 



= p[u 2 + (T - 1)6} + e PX u + e 2 (p X 2 
+pd\) 



< ; = p[u 3 + (T- l)6u] + e(p X u 2 + p(T - 1)S X ) 



,(4) 



+e 2 {px 2 u + pXSu) + e 3 {px + pX8x) 

= p[u 4 + (T- l)Su 2 + (T - l) 2 ^ 2 ] 
+e(p X u 3 + p(T - l)S X u) 
+e 2 ( PX 2 u 2 + p\8u 2 + p(T - 1)6 X 2 ) 
+e 3 (px 3 u + pXSxu) + e 4 (p x 4 + pXSx 2 ) 

where the shorthand notations of Grad [15, 74] for fully symmetric tensors are adopted. A possible set of on-site 
space-filling lattice velocities can be found in figure 12 and fully detailed in [15, 26, 30]. If one gives up the requests 
to have lattice velocities only on grid points and allows also for out of lattice discretized velocity sets, the number of 
vectors needed to recover isotropy for moments up to order eight can be reduced [30] . We next introduce [37] a small 
separation of scale parameter e and consider the expansion in e for the distribution function 



ft - h + e -ti + e Jl 



e 3 tf ' + , 4 // 4) 



and the rescaling of the time-space derivatives 

d t ^ed lt + e 2 d 2 t + 0{e 3 ); ft -> ed t . 
This allows to rewrite the streaming term in the lattice Boltzmann equation as 



(35) 
(36) 



fi(x + ciAt, t + At)- ftix, t) = cAi + e 2 A 2 
where for our purposes it is enough to consider terms up to A 2 



e 3 A 3 



p(0) 



(0). 



Ax = [ditfr + c}d lt dif^)At 

A 2 = (d 2t ft 0) + d lt f[ l) + cjdtf^At + U44didjfi 0) 



+c]d l d lt ff ) + c\d t d lt fl 0) 



d lt d lt f^)At 2 . 
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If we further rescale the shifting [37] fields as 

u^u + ex; T^T + e 2 X 
the shifted equilibrium can be further seen as a power series in e 



(37) 



f l (p,u + e X ,T + e 2 X)=f l 



(O) 



ef + ^ (2) +e 3 f 



e4/ (4) 



with 



m - P H l 



puHf' + \p\u 2 + (T - l)*]ft 



(2) 



(3) 



\2s2i-t/( 4 ) 






II 



±p[u 3 + (T - l)0u]« 

^[u 4 + (T - l)<5u 2 + (T - 1) 2 <5 2 ]^ 



,(!) 



,(2) 



(3) 



■±(pXti a +p(T-l)*x)H 
i(pXU 3 + p(T-l)5 X u)H 



(4) 



V = 5(/°X 2 + P*A)«r + \(px 2 u + P \Su)H 

-Upx 2 u 2 + p\Su 2 + p {t - 1)8 x 2 )n\ 



,(3) 



(4) 



IT 



Kpx 3 + pXd X )H^ + ±(px 3 u + p\Sxu)n ( ^ 

V 1 ^ - Upx' + P^x 2 )n\ 



// 4) 



4 ' 24* 



where, upon dimensional considerations, we have requested that when the forcing rescales as e, the temperature shifting 
term is rescaling like e 2 (see also [37] for a more detailed discussion). Using the Taylor expansion of fi(x + c; At, t+At), 
we can impose the consistency in (32) order by order in e: 



O(e ) : /, (0) 



f(0) 



i -ft 
O{^):dufl 0) +4difi 0) = -7(ti 



if) 



(38) 



J 3Jl 



(2) 



= -Hfr-fi 

Taking the momenta at the zeroth order in e we can find some constraints for the higher terms in the expansion in of 
the distribution function. Since we know that / ; = /j , it follows from the definition of macroscopic fields that 



£/," 



(") _ 



5>?jf> = o £<?//»> =o n ^ 1 - 



1. Zeroth order 



At the zeroth order in e we can find some constraints for the higher terms in the expansion of the distribution 
function. We know that 

,(0) _ j(0) 

It follows that, since we define our macroscopic variables as 



£/«; put = Y^fi c h T =D^2ft 



\ci -u\ , 



we immediately recover that 



£/r -E c w r -£i c < - u \ 2 n n) = ° »* L 
i i i 



(39) 



14 

The last equation leads to (we take the convention that double indexes are summed upon) 

% £( C W + UiU 3 ~ u%ci i _ u i cl i)fi = ° ra > 1 
I 

that, combined with the constraints for the momentum Q[\ c ;/; = 0), is equivalent to 

£c?/< B >=0 ->!• (4°) 

l 

2. First order 

We first evaluate and also remind the values of some useful quantities that can be easily obtained knowing the 
relation between Hermitc polynomials and the velocity set (33,34) and also the constraints coming from (39,40): 

£cjif =0; ^ff^PXi 



Yl $4fl = P U i U J + P TS iJ 






\ E c ^// 0) = ( W + t pT ) Ui + pTu - 



2*-^ ' '■" V 2 2 

l v 

With this, using the momenta of 0(e) in (38), we can easily arrive to the following set of equations 

' d lt p + di(pUi) = 

dit(pui) + dj(puiUj + pTSij) = 2& = gi (41) 

duJC + dj [JCuj + pTuj] = ^pXiUt = pg^ 
where we have introduced the total energy of the system: 

/c= G pu2+ f pT 

and where we have recovered the Euler equations for a forced fluid with the choice 

X = rg. (42) 

The last equation can also be written as an equation for the temperature (using the momentum equation) in the 
following form 

(du + u 3 dj)T + —T(d iUi ) = 0; c v = — . (43) 

C v L 

3. Second Order 
Using the second of (38) and the constraints found at the first order it is easy to derive: 

£ 4(d ltf w + cf5fc/ (o) } = _ i ^ cK/ (d _ /(D) = pgi _ (44) 

i i 
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Furthermore, let us write other useful quantities that can be derived from the explicit expression of the expansion of 
the equilibrium distribution, /;, and from the hydrodynamical constraints on the distribution /; reported in (39) and 
(40): 



E clcjcffj; ' = (pmUjUk + pT(6ijUk + S ik Uj + S jk Ui)) 



(45) 



E^-°; E c ^-° 



(46) 



V2i 7(i) P^'X-i , , t , D P T Xi 

?}^ c i c ifi = ~^r~ + Ui p x i u j + p TXi + — 3 — 



(47) 



9 X c ? c W// 0) = ^P^iUjU 2 - ^S ti u 2 - 



PT, 



2pTu iUj + ^DpTuiUj + ( y + 1 



(48) 






(49) 



We next proceed to evaluate some expressions in terms of the known results obtained at the previous order. In 
particular, for the momentum equation, we will have to evaluate the term: 



9u f X) c H// 0) ) = d u {pu lU] + pTSij). 



If we use the results obtained at order 0(e) in (41) we obtain 



ditipUtUj + pTSij) = -dkipuiUjiik) - Ujdi(pT) - 
Uidj(pT) + pu J g l + pmgj + 8ijpd lt T + S^Tdup 



Next, for the momentum equation, we also have to consider 



iJ„kA0 ) 



du(puiUj + pTSij) + dklpuiiijUk + 
pT(5ijUk + SikUj + SjkUi)] 



that can be simplified (with results of the previous order) as 



d u (pUiUj + pTSij) + 
+d k [pUiUjU k + pT(5 i:j u k + 6 ik Uj + Sj k Ui)] = 

pTdiUj + pTdjU % + pu, J g l + pmgj - Sij — {d k u k ). 

Cy 



For the energy equation we will have to consider 



(50) 



(51) 



9u Et c ^ (0) =9 « 



-pu 2 + —pT ) u % + pTm 



that, again, can be evaluated using the results at previous order as 



du 



Finally, we have to consider 



that gives 



f 2P u2 + —P T \ Ui + pTm 


{2 pu2 + ~2 pT ) 9i ~ dj 


ii i 



p{gkU k )u t + pTg t 



. 1 2 D 
muj [ -pu + —pT 



-2d k {pTu lUk ) - d t I -pTu 2 \ + pT Uj d t u 3 

— + 1 J Tdi(pT) pTui(d k u k ) + puidjUj. 



<±A t(°) I _l a. ir <±AJ t(P) 



ftilEf^J+^IEf^/. 



1 , D 

-pu + — pT ) m + pTui 



pUiUjU H (Jyu + 2pTuiUj 



-DpTuiUj 



D 



1 pT 2 <5, 



feETtf MEW 
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(52) 



+p{9kU k )u l + pTg t 4 ( -pu 2 + —pT ) g t -+ 



^ + l)pT dl T + pT(u ldj u 3+ u A u^ 
pTui(d k u k ). 



(53) 



We are now ready to write down the equations at this order using results in (45)- (48) and (50)- (53) 

'd 2t p + l/2di(p gi At) = 

d2t(pui) + djirpgiUj + Tpgjm) + d lt (^M) = 



Af 
2 
.PT 



dj{pTdiU + pTdjU, + pu j g l + p mgj 



EhtlC + d k [JCTg k + rpu k g k + pTg k ] + d u ( £S 2 iii Ai) - 
(r - M) d l [p{g k u k )u l + P T 9l + K gi + (f + l) pTd{T 
+pTu 2 djUj + pTujdtUj - j-pTui(d k u k )] = 

{=H\pr 2 g 2 + \D P x). 

Summing up all orders, we note that wc can freely add at order 0(e 2 ) all the gradients of terms 0(g 2 ) and also double 
gradients of terms O(g) because they would be 0(e 3 ). Also, defining the hydrodynamic velocity as u. 



(H) 



Ui 



gjAj 
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we reconstruct the following equation s: 

'd tP + d i (pu\ H) ) = o 



d t (pui H) ) + diifm^uf*) = -di(pT) + 9i 

) d 3 \pTd iU ( H) + pTdjU^ - S^dku^ 



An 



d t lCW +d j \lCWuf I) + pTuf 
}- T (pr 2 g 2 + DpX) + (t-%) di[(±D + l) P Td t T 

v -ou {H) (diu {t 



pTu\ H) d lU f ] + P Tu ( " ) d l u (H) 



= P9kUk+ 



with 



In order to recover the correct thermohydrodynamical evolution we need to obtain the correct forcing in the equation 
for the total energy in terms of the hydrodynamical velocity fields, i.e. 



P9ku k + — {pr 2 g 2 



DpX) 



(H) 

P9kU k 



P9k u k 



Atg h 



that leads to 



A 



T(At - t) 9 2 
D 



(54) 



In conclusions, expressing everything in terms of the hydrodynamical fields, it is easy to realize that the final expression 
(54) coincides with the one given in the body of the article (9). Notice that up to now we have used a single-time 
relaxation LBM, as given by (32). Therefore, the final Fourier-Navier- Stokes equations are constrained to describe 
fluids with unit Prandtl numbers, Pr = v/{k/c p ) = 1. It is possible to generalize the approach by using a multi- 
relaxation time version of the same algorithm [29] . Even though, in the latter case, there exists a small mismatch in 
the viscous dissipation term appearing in the energy balance. 



VII. APPENDIX B 



In this appendix we detail the technical steps leading to the desired hydrodynamical boundary conditions for the 
physical systems analyzed in the paper, i.e. an ideal gas under the effect of gravity g = (0, — g) acting along the 
negative z direction (i.e. g is positive). Similar ideas can be applied to the case of a generic volume or internal 
force acting also in the stream-wise x direction. For the sake of concretcness we explicitly report the case of the 
lower boundary condition with the upper boundary condition being a straightforward generalization. Let us call the 
post streaming populations /* while keeping _/y ' p to identify the pre streaming populations. Moreover, all the 
populations will also undergo collisions and therefore there will be a net gain of momentum so that the hydrodynamic 
fields will be the average of pre and post collisions. For a given computational boundary, there are 3 layers of 
points labeled by x* from now on (see also figure 13), where some unknown populations have to be set soon after 
the streaming step. We use the freedom to set these populations in such a way that the measured hydrodynamic 
quantities such as the stream- wise {u x ) and vertical (u z ) velocity and also the temperature (T^) are fixed to 
some given boundary conditions on those lattice layers. The conditions to be fulfilled up to the second order in the 
Chapman-Enskog expansion (see also previous appendix) are 



,W 



GO 



p(x 



^rE/fOO 



c ) 



(55) 



iW 



oo 



p(x 



^£/roo 



c i 



At 



(56) 
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E T m( x *) + l((uWf + (vWf)(x*) = 



2p{x*) ^ 



f!{x*)cl (57) 



In the following we show how to determine the unknown populations on the first three layers (those coming -after 
streaming- from node outside the domain) in order to set the vertical velocity to zero on layer 3, with any temperature 
and stream- wise velocities: 

~ui H) (z = 3) = Q 

Ux (z = 3) = U3 

TW(z = 3)=T 3 . 
Similarly we can fix any desired profile for temperature and velocity on layers 1,2: 



a 



(H \z = 2)=v 2 - u{ H \z = l) = v 1 

(H),_ n _„, . (H) 



u K x >(z = 2)=u 2 - 1 Mi > (z = l) = u 1 

TW{z = 2)=T 2 ; TW(z = 1) = T 1 . 

We will define only the case of homogeneous boundary conditions along the stream- wise component but the method 
is general and can deal also non-homogeneous cases. Imposing a given set of boundary conditions means defining 
the set of unknown outgoing populations in the first three layers in terms of the set of in-going and outgoing known 
populations such that mass is conserved and the hydrodynamical fields defined above are the wanted ones. 
In this way, if the computational boundary extends from the mesh point z — 1 up to z = L z , the real physical domain 
is between mesh points z = 3 and z = L z — 2, i.e. it is in these points that we exactly verify the condition of no-slip, 
no normal velocity and given temperature for the hydrodynamical fields on the solid walls. Fields at points z = 1, 2 
and z = L z — 1, L z — 2 may be used to better stabilize the algorithm close to the boundaries. All details refer to the 
37 speed model D2Q37. 

Layer 1 

As evident from figure 13 we have to determine some 'outer' post streaming populations (I = 2, 10, 18...) whereas 
other post streaming populations (/ = 4, 12, 20...) are known. To keep a compact notation, let us also introduce the 
subsets I 1 - 1 ^, U^ and Iq which arc identified by the following conditions 

lW={c h cf <0}; U^={c h c?>0} 

I^ = { Cll ci<i)}. 
We choose to define the 'outer' populations in the layer 1 as 

/i (1 '" = — JL -m^ ) leffCD (58) 

with N a constant and 4>l a suitable population that we choose in the form 

^^l + C-pW + ^l (59) 

where p x ,p z and E^ l > are unknown at this level and must be chosen in such a way that the hydrodynamical 
temperature and momentum exactly reproduce the desired values on this layer, Ti,ui,V\. Also, mass conservation 
should be fulfilled. This latter condition is naturally imposed by setting 

N= ]T f^*' pre \ 
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The requirement that T\,u\, Vx are exactly reproduced leads to the following system of equations 



^ = ^Ejrcf-#ff (60) 



where we have defined the post streaming mass as 

In the ^ ; fi of system (60) we have known populations coming from the bulk but also 'outer' populations to be 

determined with (58) and (59). The resulting system is therefore an algebraic system for p x ,p z ' and E^\ We have 
solved the system whose final solution is 

(i) -c 3 d 2 + c 2 d 3 

p\ = 

~a 3 c 2 + a 2 c 3 

(i) _ 0203^1 — 0201^3 — a 3 c 2 di — c 3 a\d 2 + c\a 3 d 2 + a\c 2 d 3 

Px 



where 



with 



and 



bi(a 3 c 2 - a 2 c 3 ) 

(1) _ -a 3 d 2 + a 2 d 3 
a 3 c 2 - a 2 C3 



oi = 26(p x - O x )r- bj. = -AONr 2 
c 1 =A7{p x -O x )r 2 ; dx = 15(p x - O x ) 
a 2 = 26{p z - O z )r - 5ANr 2 ; c 2 = 47{p z - O z )r 2 - 9lNr 3 
d 2 = 15(p 2 - O z ) - 26iVr; a 3 = 26{E - O e )r - 9lNr 3 



367 
c 3 = 47(£ - O e )r 2 - —Nr 4 ; d 3 = 15(JS - O e ) - A7Nr 2 



p{ 1] = M pUl ; p[ 1] = M p vx + ^M p gAt 



E = T 1 M p + ^-((p x V) 2 + m 2 ) 



O x = £ eff^; O z =J2 ^// M) 



In the above r is the lattice constant whose value for the D2Q37 model is r ~ 1.1969 [29] 
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Layer 2 

Situation goes similarly with respect to the previous layer. We new have to define the subsets I^ 2 \ U^ and Iq as 

I^ ={c u cf <-r}; U^ = {c h cf>r} 

4 2) = {Q,cf<r}. 
We then identify some coarse grained quantities as 

N= ^/P'* 3 ; M P = N+ J2f^ ] 
leim leI m 

ad define some local momentum and energy fields 

p x 2) = M p u 2 

pW = M p v 2 + -M p gAt 



We next define 



= -!■'. 


>Mp -\ 


2Mj 


i 


+ u 


?rj ;• 


E 


r x f('' 
C l J I 




o,= 


E 


c Ul 


'€tf J 








^ 2) 






o e = 


E 


1 r 2 f C 


!,*) 





oi = 19 (p x - O x )r; &i = -12iVr 2 
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ci = yfe -O x )r 2 ; dl=8{p x -O x ) 



a 2 = 19(p, - 0,)r - 477Vr 2 ; c 2 = y (p, - O z )r 2 - ^Wr 3 
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^2 = 8(p 2 - 2 ) - 197Vr; a 3 = 19 (E - O e )r - —Nr 3 



C3 = f (£ - O e )r 2 - ^iV r 4. ^ = 8( £ _ 0e) ^ Nr 2 



In terms of these constants and parameters we can set 

( 2 ) _ -c 3 d 2 + c 2 d 3 



Pz 

— a 3 c 2 + a 2 C3 



(2) _ «2C3tii - 02^1^3 - 03^2^1 - c 3 aid 2 + cia 3 d 2 + a\C 2 d 3 

Px 



bi(a 3 c 2 - a 2 c 3 ) 
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-a 3 d 2 + a 2 d 3 



03 C2 - CI2C3 
construct suitable populations 



-1 + c.pW + i^W 

and define the outer populations in the layer 2 as 



»l = J-t-ci-p- • -i--q. 



that is enough to set the hydrodynamic velocity to u 2 and v 2 while keeping the hydrodynamic temperature to T%. 

Layer 3 

As also evident from the figure 13, only 3 populations are unknown on the third layer (they are populations 
/ = 24,25,18). In this way we do not have enough freedom to choose the desired hydrodynamic velocities and 
temperature. It is anyhow possible to require a zero vertical hydrodynamic velocity (113 = 0) with a generic stream- 
wise hydrodynamic velocity and temperature (1*3^3). Again, let us introduce the following sets 

U^={ Cll cf >2r}; /f = {c h cf < 2r}. 
The boundary condition for the unknown populations is set as 

^-1 + cfp^ + lcfE^ 

and we choose p x and E^ 3 ' to set the desired hydrodynamical stream- wise velocity (113) and temperatur e (T3) while 
keeping the vertical hydrodynamical velocity to zero. The resulting algebraic system is solved with the solution 

(3) _ _d2^ (3) 6id 2 -rfi6 2 

CJ ■ - ' Px 



62' e * ai6 : 



with 



2 



2. ,29,, „,, 



where 



ai = -2Nr 2 - b l ^—{p x ~O x )r 



d, = 3(p x - O x ); b 2 = ^-(E-O e )r 2 -^-Nr 4 



29 
d 2 = 3(E - O e ) - —Nr 2 



ze/< 3 > ;e/< 3) 



n(3) 



pW = M p u 3 



^ (3) - r 3M p + -i-((pi 3 )) 2 + (p( 3) ) 2 ) 



2M p 
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pW=3Nr+ Vcf/.^ 






M P = iV + J2 fi 



'e-A 



(3) 




4i 1 A 2 A 

TV = + -gA£ 

3r 2 3r y 



* = E c ^ (M ; ^= E/f^-if^ (6i) 



ze/i ; zei. 



(0) 



This whole algorithm for layer 3 now is ensuring a zero vertical hydrodynamical velocity and arbitrary U3 and T3. 
Still, mass conservation is not fulfilled and to do that we need to redefine the rest population appearing in (61) as 



f (S,*) =f (S,*,pre)_ N+ J2 /< 
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run A 
run B 



/i-t Lix 



Lad 



T u T d N, 



0.05 800 2400 4 x lO 3 0.001 5 x 10~ 5 0.95 1.05 
0.4 1664 4400 1.6 x 10 4 0.1 1 x 10~ 4 0.6 1.4 



conf T 

50 1.8 x 10 4 
35 6.5 x 10 3 



TABLE I: Parameters for the two sets of Rayleigh- Taylor run. Atwood number, At = (Td — T u )/(Td + T u ); Adiabatic Length, 
Lad — (Td — T u )c p /g (c p = 2); viscosity v; gravity g; temperature in the upper half region, T u ; temperature in the lower half 
region, Td\ number of separate Rayleigh- Taylor run N co „f, normalization time, f = \fL x /(g At) (not to be confused with the 
relaxation time of the lattice Boltzmann model (2)). Given the parameters here used, the typical resolution obtained is good 
enough to get an agreement better than a few per cent on the global exact balance between kinetic energy growth and the sum 
of dissipation plus buyoancy force. 



27 



L z /2 



-L z /2 




L z /2 




FIG. 1: Upper panel: Rayleigh-Benard geometry and set-up of the initial configuration given by eq. (14); two cases with 
m — +0.5 and m = —0.9. On the horizontal axis we show the mean temperature and density profiles as a function of the 
z- height (plotted on the vertical axis). The bold and tiny solid lines represent the temperature and density profiles respectively. 
Lower panel: Rayleigh- Taylor initial configuration given by eq. (17). Bold and tiny lines as in the upper panel. 
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FIG. 2: Critical Rayleigh number, estimated at the center of the cell, Ra c , at changing the polytropic index, m, the scale 
height, Z, and the numerical resolutions. For the smallest resolution, L x x L~ = 232 x 116, the plotted values corresponds to 
(a) Z = 0.22, AT = 0.2, T d = 0.9, m = -0.942; (b) Z = 0.5, AT = 0.4, T d = 0.8, m = -0.971; (c) Z = 0.86, AT = 0.6, T d = 
0.7, m = -0.9806; (d) Z = 1.33, AT = 0.8, T d = 0.6, m = -0.9855; (e) Z = 2.0, AT = 1.2, T d = 0.6, m = -0.990. Theoretical 
values are obtained solving the linearized equations as described in [34]. Inset: time evolution of the total kinetic energy (in 
arbitrary units) for Rayleigh numbers lower and higher than the critical one for the parameter case (c). The unit of time 
corresponds to 10000 LB integration steps. 
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FIG. 3: Comparison between one-dimensional vertical cut of the stationary temperature and velocity profiles after transition 
to a convective two-rolls configuration. Up: T(xo, z) at xo such that xo = 0.69L X . Circles correspond to the Lattice Boltzmann 
algorithm (LBM); solid line corresponds to a finite difference calculations (FDM) [55, 56]. Down: the same of above plot but 
for the stream- wise velocity, u x (xq, z). In the insets we show a grey-scale coded representations of the convective stationary 
rolls in the whole two-dimensional domain (up: temperature; down: stream- wise velocity). 
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FIG. 4: Spatial configuration for a typical RT run with L x x L z — 800 x 2400, T u = 0.95, T d = 1.05 at time t — 4f (run A in 
table I). 
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FIG. 5: Temporal evolution of the mean temperature profile, {T(x,z,t)) x at changing time, t = nSt, with St ~ l-5r a d 
n = 1,2,..., 7. Notice that the profile approaches more and more the linear behaviour dictated by the adiabatic gradient, 
{T(x,z,t)) x = (T u + Td)/2 — zg/cp. Time is adimensionalized by using a reference time based on adiabatic quantities, r a d = 

VW(<?AT/T U ) 
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FIG. 6: Evolutions of the mixing layer, L m i(t) versus time with two different adiabatic lengths: (a): L a d — 800, g — 2.5 10 4 
(triangles); (b) L a d — 400, g = 5 10~ 4 (circles); Both cases have At — 0.05, v — 0.001 and k, — 0.002. Solid bold lines correspond 
to the theoretical prediction (23) with or- ' = 0.05. Continuous line correspond to the evolution of the instantaneous Rayleigh 
number (20) calculated for case (a), scale on the right y-axis. 
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FIG. 7: Temperature, (T(x,z;t)) x , density, {p(x , z; £)) x , and pressure, (p(x,z;t)) x , instantaneous mean profiles at different 
time during the RT evolution. Left column: At — 0.05, times t = 3f , 6f , 7f (run A, table I); Right column: At = 0.4, times 
t = 3f, 4.5f,6f (run B, table I). Initial hydrostatic profiles are depicted by solid bold lines. Notice the asymmetry for the 
mixing layer growth in the latter case. Notice also the appearance for high Atwood of hydrodynamical pressure fluctuations 
superposed to the hydrostatic pressure profiles. Both effects are absent in the small Atwood case. 
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FIG. 8: Growth of H„ : b(t), run-by-run, for the two Atwood numbers of run A and run B. The mixing length width is normalized 
by the total box width, L z . Notice the evident asymmetry between spikes and bubbles for the high Atwood case (run B). 
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FIG. 9: Run A. At — 0.05. Analysis of the asymptotic growth rate for spikes, a s . Bottom panel: mixing length evolution 
normalized by t 2 for two different set of initial width, W = e/uj — 2 (circles) and W = 8 (triangles) , where e is the intensity 
of the initial perturbation and lo is the width of the regularizing tanh initial profile. Data refers to N con f = 50 for both 
cases. Notice the long relaxation time before the two evolution forgets the initial conditions. This is due to the presence of the 
prefactor proportional to L m i(to) in the sub-leading linear term of eq. (22). Inset: the mean value of the data shown in the 
body but in log coordinates -same symbols. 

Upper panel: mean value of the instantaneous growth rate of spikes extracted from (30) for the two initial set up with W = 2, 8. 
Average is performed over N con f = 50 separate Rayleigh Taylor evolution for the two cases. Error bars are estimated out of 
root mean square fluctuations. Notice the more extended range where the two set-up superpose and the extended time interval 
where a a stays constant (notice the different y-scale between lower and upper panels). Results for bubbles evolution are 
similar and not shown. Both cases are summarized in figure 10 
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10: Top: Run A, At = 0.05; histograms of a s b (multiplied by 10 3 for the sake of clarity) as extracted from (30), at fixed 



initial width W = 2 (left) and W = 8 (right). The fit is done over 50 and 35 different configurations respectively. In order to 
test dependency on the fitting window we have summed results from two different ranges, t € [1.5f : 4.5f] and t € [2.2f : 4f] in 
both cases the maximum time is such that the front didn't reach more than 80% of the total vertical extension of the physical 
domain. Bottom: Run B, At = 0.4. Results from two fitting ranges t £ [2.3f : 5.4f] and t € [3f : 4.5f]. Notice the asymmetry 
developing for At = 0.4, with spikes traveling faster. An estimate of the mean value for the growth rate in the two cases gives: 
ai H) = (10 ±2) 10~ 3 and a l b H) = (9.5 ±2) 10" 3 at At = 0.05, while a { s H) = (14 ±4) 10" 3 and af = (9 ±5) 10~ 3 for At = 0.4. 
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FIG. 11: Bottom: 2nd order moment of temperature fluctuations (see Eqn. (31)) as a function of the height at times 
t = (2.2, 3.3, 4.4)f . The z-height has been normalized with the total cell extension L z . Top: Flatness Q' 4 '(z, £)/(Q' 2 '(z, t)) 2 at 
the same three instants of time as in the bottom panel. The z-height has been normalized with L m i(t) in order to show the 
self-similarity of the mixing process (the three curves collapse onto each other by rescaling). Parameters refer to run A in table 
1. 
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FIG. 12: Scheme for the D2Q37 model used for the simulation of thermohydrdoynamics. The 'lattice constant' is r ~ 1.1969 
as reported in [29]. The velocity set is such that every projection of the velocity is an integer multiple of r which is chosen 
to enforce the unitarity of Hermite polynomials (33,34) up to the fourth order. The relationship between real and velocity 
lattice is set by Ax = rAt with Aa; and At space and time discretizations. Based on the Hermite-Gauss quadrature procedure 
[15, 26, 29], the D2Q37 can be regarded as the minimal on grid square lattice giving with accurate Hermite polynomials up to 
the fourth order. This quadrature ensures that the Navier-Stokes thermodynamics is recovered with full Galilean invariance. 
Lattice D2V37 firstly appeared and was shown to be minimal for 2d fourth order models in Reference [26], where the authors 
formally showed the equivalence between the condition of norm preservation and the preservation of the orthogonality property, 
in constructing these sets of lattice vectors. 
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FIG. 13: Scheme for the lower boundary layer for the simulation of thermal flows under the effect of gravity. The relationship 
between real and velocity lattice is set by Ax — rAt with Ax and At space and time discretizations, and r the lattice constant 
whose value is r ~ 1.1969. the locations at r and 2r indicated in this figure correspond to the locations 2 = 2 and 2 = 3 
discussed in the text. 



